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Abstract: We study the perturbative behavior of the Yang- Mills gradient flow in the 
Schrodinger Functional, both in the continuum and on the lattice. The energy density of 
00 , the flow field is used to define a running coupling at a scale given by the size of the finite 

volume box. From our perturbative computation we estimate the size of cutoff effects of 
■^-j- \ this coupling to leading order in perturbation theory. On a set of Nf = 2 gauge field 

ensembles in a physical volume of L ~ 0.4 fm we finally demonstrate the suitability of 
the coupling for a precise continuum limit due to modest cutoff effects and high statistical 
precision. 
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1 Introduction 

Finite-volume renormalization schemes have now a long history in lattice field theory 
(see [1-3] or the pedagogical reviews [4, 5]). Asymptotic freedom tells us that at small 
distances QCD is well described by perturbation theory, while at large scales QCD is a 
strongly interacting theory. Instead of trying to accommodate these two scales in a single 
lattice simulation, the idea of finite-size scaling exploits the size of a finite volume world as 
renormalization scale. A single lattice simulation covers only a range of scales, but one can 
match different lattices and adopt a recursive procedure to cover a large range of scales. 
In this way one can connect the perturbative and non-perturbative regimes of QCD. 
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Beside a successful application of the finite-size scaling technique in case of pure Yang- 
Mills theory [1-3], the running coupling [6-10] and quark mass [11-13] have been computed 
non-perturbatively in QCD with different flavour content by ALPHA and other collabo- 
rations. Since the general idea of finite-size scaling is a very powerful tool to solve scale 
dependent renormalization problems, it is not surprising that it is broadly used also in other 
strongly interacting theories, even in effective theories such as HQET [14-16]. There has 
been a growing interest in other than QCD strongly interacting gauge theories, especially 
in connection with electroweak symmetry breaking and quasi-conformal behavior (see for 
example [17] and references therein). Finite-size scaling techniques are also a powerful tool 
to study these systems. 

Basically there are two things that are needed to perform the previously sketched 
program. First one needs to define exactly what is meant by a finite-volume scheme, 
i.e., one has to specify the boundary conditions of the fields. Second, one needs a non- 
perturbative definition of the coupling. In principle there are many valid possibilities, but 
practical considerations have to be taken into account. Good options should allow for an 
easy evaluation of the coupling constant both in perturbation theory and in a numerical, 
non-perturbative (lattice) simulation. 

The rest of this section is mainly dedicated to explain why we choose the Schrddinger 
functional (SF) scheme [2] as our finite-volume setup and the Wilson flow for a non- 
perturbative definition of the coupling [18]. To simplify the following discussion we will 
argue about a pure SU(N) gauge theory in 4-dimensional Euclidean space-time. 

In the Schrodinger functional [2, 19] one embeds the fields in a finite volume box of 
dimensions L 3 x T. Gauge fields in the SF are periodic in the three spatial directions and 
have Dirichlet boundary conditions in time direction (i.e. one fixes the value of the gauge 
fields at xo = 0, T) . The value of the gauge fields at the time boundaries are called boundary 
fields. One can interpret the partition function of the theory as the probability amplitude 
of the gauge field to propagate from the boundary value at xq = to the boundary value 
at xo = T. Such a setup has nice properties in perturbation theory. In particular with a 
smart choice of the boundary fields one can guarantee that there is a unique gauge field 
configuration (up to gauge transformations) that is a global minimum of the action. This 
avoids some difficulties with perturbation theory [20-22]. The reader interested in this 
issue would appreciate the original literature, as well as the nice discussion in [23]. 

Recently, the gradient flow has been used in different contexts [24-26], but it is the 
proposal made in [18] to define a renormalized coupling through the gradient flow in non- 
abelian gauge theories what inspires this work. The gradient flow defines a family of 
gauge fields parametrized by a continuous flow time t. The flow equation brings the gauge 
field towards the minimum of the Yang-Mills action, and therefore represents a smoothing 
process. The key point is that correlation functions of the smoothed gauge field defined at 
t > are automatically finite [27]. One can use the expectation value of the energy density, 



where G^ v {t) is the field strength of the gauge field at flow time t, to give a non-perturbative 



{E{t)) = -{G, u {t)G, u {t)), 
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definition of the gauge coupiing. This idea was applied to set the scale in lattice simula- 
tions [18, 28], to tune anisotropic lattices [29] and more recently in a similar context of this 
work (finite-size scaling, but using a box with periodic boundary conditions) to compute 
the step scaling function in SU(3) with four fermions [30]. 

In this paper we investigate the perturbative behavior of the Wilson flow in the 
Schrodinger functional. This motivates us to propose a gradient flow coupling 

f GF (L) = N-h 2 {E(t)) = g 2 MS + O^s), (1.2) 

with a normalization factor M to be determined later, valid for an arbitrary SU(N) gauge 
field coupled (or not) to fermions. The coupling depends only on one scale, the size of the 
finite volume box, and therefore can be used for a finite-size scaling procedure in the same 
way as the traditional SF coupling. 

The paper is organized as follows: in the next section we investigate the perturbative 
behavior of (E(t)) in the SF, both in the continuum and on the lattice. Section 3 uses this 
information to define the gradient flow coupling in the SF, and to discuss some practical 
issues: cutoff effects, boundary fields and fermions. In section 4 we investigate this coupling 
numerically on a set of lattices in a physical volume of L ~ 0.4 fm and finally conclude in 
section 5. Details needed for the computation have been summarized in form of appendices: 
a summary with some useful notation A, heat kernels B, propagators in the SF C and finally 
some practical details on how to integrate the Wilson flow in numerical simulations D. 



2 Perturbative behavior of the Wilson flow in the SF 

We would like to start this section by recalling the original proposal of using the Wilson 
flow and the energy density as a definition for a coupling in gauge theories [18]. Later it 
will become clear what role the SF setup plays. 

2.1 Generalities 

By considering the gauge fields to be functions of an extra flow time t, not to be confused 
with Euclidean time, denoted xo, the Wilson flow is defined by the non- linear equation 

dB ^ f) =D„G^(x,t), B^(x,0) = A ll (x), (2.1) 

where 

= d^B v - 8 U B^ + [B^ B v \ . (2.2) 

Due to D u G Ufl ~ 'sb " — " S^ug^ fields along the flow become smoother, eventually 

reaching a local minimum of the Yang Mills action: the flow smooths the fields over a region 
of radius y/St. The somewhat surprising result of [18, 27] is that correlation functions made 
of this smooth field have a well-defined continuum limit. In particular the energy density 
in SU(N) Yang-Mills theory in infinite volume has the perturbative behavior 

(E(t)) = \{G, U G, V ) = 3( ^ 28 ~ ffi S (1 + 4 S + 0{tus)) , (2-3) 
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where c\ is a numerical constant and g MS is the renormalized coupling in the MS scheme. 
Therefore one can define a running coupling constant a(p) at scale \i = \j\f%i from 

t 2 (E(t)) = 3( ^ 2 ~ [1 + <Wm) + 0(« 2 )] ■ (2-4) 

These expressions are valid in infinite volume. What about the Schrodinger Functional? 
The computation is completely analogous, but we have to impose the correct boundary 
conditions to the gauge fields. As we have mentioned in the SF gauge fields are embedded 
in a box of dimensions L 3 x T. They are periodic in the three spatial directions and have 
either Dirichlet or Neumann boundary conditions in the time direction. We are going to 
work exclusively with zero boundary fields, which means 

Bp{x + kL,t) = Bn(x,t), (2.5) 
B k (x,t)\ xo=0 ,T = 0, (2.6) 
d Bo(x,t)\ Xo=0 ,T = 0. (2.7) 

The flow equation (2.1) has to be solved maintaining these boundary conditions at all 
flow times t. To apply the idea of finite-size scaling, as has previously been done in [23] in 
a periodic box, one simply has to run the renormalization scale with the size of the finite 
volume box given by L via 

« = Wt = k- < 2 - 8 > 

Here c is a dimensionless constant that represents the fraction of the smoothing range over 
the total size of the box. In this way the flow coupling will not depend on any scale other 
than L. The renormalization scheme will depend on the values of c, p = T/L and 1 xq 

g 2 GF (L)=N~ 1 (c,p,x /T)t 2 (E(t,x )) , (2.9) 

t=c / L / /8 

where Af (c, p, xq) will be computed in the next section in order to ensure 

5gf = 3o + 0(3o)- (2-10) 

2.2 Continuum 

Our computation follows the lines of [27]. By using a flow time dependent gauge transfor- 
mation, the Wilson flow equation can be casted into the form 

^■ = D v G„ li + D fl d v B v . (2.11) 

After rescaling the gauge potential with the bare coupling — > goA^, the flow becomes a 
function of the coupling 

B fl (x,t)=J2 B i^n(x,t)gS. (2.12) 



1 Note that in the SF the boundary conditions break the invariance under time translations. Therefore 
(E(t,xo)) will depend explicitly on xq. 
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Inserting this expression in the flow equation, we find that to leading order in go the flow 
equation is just the heat equation: 



dt 



d 2 u B^(x,t), (2.13) 



i.e., to leading order the Wilson flow is the heat flow. The heat equation can be solved by 
using the heat kernels 

r-L r-T 3 

B M (x,x ,t)= / dV / dx' Q T\K p (x i ,x' i ,t)K D (x ,x Q ,t)A k (x , ,x ), (2.14) 
Jo Jo i=l 

r-L r-T 3 

B 0)1 (x,x ,t) = / d 3 x' / dx' r[K p (x i ,^,t)^ 7V (xo,x , ,t)^o(x , ,4). (2.15) 
Jo Jo l=1 

Since the boundary conditions of the field B^\{x., xq, t) are inherited from the boundary 
conditions of the heat kernels, we have to choose them with correct boundary conditions. 
In particular 

K p (x,x',t) = lj2e-*>\^-*'\ L = ^) (2.16) 

is the one-dimensional heat kernel with periodic boundary conditions on [0, L\. K D (x, x', t) 
and K N {x, x' , t) have Dirichlet and Neumann boundary conditions respectively on the 
interval [0, T\. See appendix B for explicit expressions and more details. 

Our observable, the energy density (E(t,xo)), has an expansion in powers of go- The 
leading contribution is given by 

So(t,x ) = fid^d^ - d^dvB^) . (2.17) 

We are going to split the computation in two parts, one involving only the spatial compo- 
nents of G^y , and the other involving the mixed time-space components of G M „ 

£ s o{t,x Q ) = ^(diB^diB^ - diB^dkB^) , (2.18) 
£o(t,x ) = ^{doB^doB^ - doB^dkB^) . (2.19) 



By expanding the gauge field 

A^) = ^E eipx VP' x °) ( 2 - 2 °) 

p 

and inserting the heat kernels with appropriate boundary conditions in (2.17) we obtain 
m^o) = - ^e^ 2+ ^h<^ T dx'ody'oK D (xo,x'o,t)K D (xo,y'o,t) 

pq J ° 

x [mi(A%(p,4>)A%Mo)) -Pi^(A a (p^d)^(q,2/o))] • (2.21) 
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The final result is obtained inserting the SF gluon propagator [31, 32], that in the Feynman 
gauge reads (for additional details see appendix C) 

(AUP, *o)4(q, yo)) = LX b 6 lk 5± £ f^lfgif) . ( 2 .22) 

1 po P 2 +(f) 



To shorten notation we use 



s P0 (x) =sin(^) , (2.23) 



c po (x) = cos(— J , Po = ^^. (2.24) 



After some algebraic work one arrives at the expression 

•A ( AT2 1 "\ „2 97/9,1 9\ ^2 



tt Q {t,x ) t g - 64 2^ e n2+ iJ^W) ^ 25 ^ 

' n,no 4p z U 



and a very similar computation leads to 



t£o(t,x J g - 12g 2^ e n 2 + i n 2 c no^o). (2.26) 

r n,no 4p^ U 

2.3 Lattice 

On the lattice one defines the Wilson flow as 

a 2 d t V^x,t) = -<7 2 {T^^(y)}^(x,t) , V^(x,0) = C/^x) , (2.27) 

where V^(x,i) are the lattice links associated with the flow field Bu(x,t), Uu(x) those 
associated with the initial gauge field A /Ji (x), and S w is the Wilson action. If f(U^(x)) is 
an arbitrary function of the link variable U^x), its Lie-algebra valued derivative <9"^ is 
defined as 

%,„/( W) = d/(e£T ^ (x)) . (2.28) 

To first order in <?o and after a flow time dependent gauge transformation, the lattice 
flow equation reads 

dtB^x {x,t) = 8 u dtB^ (x, t) (2.29) 

where d and d* are the usual forward/backward finite differences (see appendix A). Now 
we have to solve a special type of heat equation in which the Laplacian is substituted by 
a discrete version, but the flow time remains a continuous variable. The strategy is very 
similar: We find the fundamental solutions of this equation, i.e., the discrete heat kernels 
given in appendix B, and write 

L-l T 3 

B^{-K,x ,t) = Y,Y1 nK P {xi,xlt)K D > N (x Q ,x' Q ,t)A^ ,x' Q ) . (2.30) 

x'=0ij=0 i=l 
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Then we have to insert this in our lattice observable (E). We use the clover definition for 
G^y that to leading order in go reads 

G>„ = y d, t [B Vy i(x) + B Vtl {x - i>)] - y 8 V \B^{x) + B^{x - p)] + 0(gl) , (2.31) 

where = \{d^ + d*). On the lattice the gauge field is expanded as finite sum over all 
lattice momenta according to 

Mx) = ^E e %P ' X e mPk/2 MP, *o) , (2-32) 
p 

Mx) = ^3^e i P x i (p,x ), (2.33) 
p 

and the computation is completed by using the lattice gluon propagator 

(Af(p, so)ifc(q, yo)) = LH ab 5 lk 5^ 1 Y, ^f" P i V " ] ■ (2.34) 

po u 

In appendix A we give the definition of the lattice momenta p, p, p and the functions 
Sp(x),c p (x). For the spatial part of the contribution to the energy density we arrive at 

0V U; i=c 2L2/8 128/9 ^ 

P,PO 

p 2 cos 2 (api/2) - (pi cos(api/2)) 2 „ 2 

- 2 , - 2 Spo^o ( 2 -35) 

P + 

while for the mixed part we obtain 

t^(t,x ) = ^Ve-^(P 2 ^ x 

v u; i= c 2 L 2 /8 i28p ^ 

P,PO 

p 2 cos 2 (apo/4) + \p\ cos 2 (api/2) , 2 

~ 2 ~ . 2 c po x o - a/2) . (2.36) 

P 2 + Po 



2.4 Tests 

There are several tests that can be performed to check the previous computations. Let 
us first concentrate on the continuum computation. At fixed t, the infinite volume limit 
L — > oo (with p kept constant) is taken through c — > 0. For this case the continuum 
expression transforms into an integral via 

cn — > p , -k — > po , — V — > [ d 4 p, 
P P ^ J 

n,k 

and we obtain 

lim[^(,T/2) + ^(,T/2)] = / ^ Pe ^^) ^ + P + 3P8 

-MO (2 37) 

" 128^ 2 ^' d7j 
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(a) c = 0.3 and xq/cl = 6. 



(b) c = 0.4 and x /a = 4. 



Figure 1: Two examples of the perturbative prediction versus a large f3 pure gauge sim- 
ulation. The fit to a linear behavior in g 2 intercepts at zero with a precision of about 
10" 4 . 



thus recovering the infinite volume result of [27]. 

Another rather obvious check is that one should recover the continuum result from 
the lattice expression in the limit a/L — > 0. This can be easily checked by noting that the 
sums (2.35) and (2.36) are dominated by terms with small lattice momenta ap^ — > 0. 

Finally we have performed some simulations with the openQCD code [33] at small values 
of the bare coupling in a pure SU (3) gauge theory. Using a 8 3 x 7 lattice and varying 
the value of the bare coupling (/3 = 60, 120, 240, 360, 480, 600, 720, 840, 960, 1080, 1200) we 
compare the analytical lattice prediction and the numerical results after collecting 10000 
measurements of the gradient flow coupling for each value of f3. We use the clover definition 
for G^y to compute the value of 

t 2 (E(t,x ))\ t=c 2 L2/s . (2.38) 

The lattice computation of 

t 2 £ (t, s„) = t 2 xo) + 4 m (t, s )] (2.39) 
can be checked in the following way: plotting 

(E(t,x )) - £o{t,x 



L2/8 = O(g 2 ) (2.40) 



£ (t,x ) 

versus g$ one expects a linear behavior with zero intercept for all values of c and xq. A 
couple of typical cases are show in the figure 1, while table 1 shows the results of the x 2 /dof 
of the fits and the intercepts for all values of c and xo . 

All intercepts are of the order 10 -4 and compatible with zero within errors. The 
difference in £q{c, xq) for different values of c, xq or between the continuum and the lattice 
result varies between 5% and 10%. We note that this last test is highly non-trivial since 
it is done for arbitrary xq at p ^ 1 on a small lattice where finite-size effects tend to be 
larger. 
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c 



xq = 1 xq = 2 xq = 3 xq = 4 xq = 5 xq = 6 



0.3 



1.77 1.53 1.31 1.68 2.10 1.70 
5(7) -2(5) -6(5) -8(5) -4(6) 4(6) 



X 2 /ndof 
Intercept x 10 4 



0.4 



0(1) 0(7) -4(5) -7(6) -3(8) 5(8) 



1.37 0.98 0.89 1.35 1.34 0.91 



X 2 /ndof 
Intercept x 10 4 



Table 1: Parameters of the fits to the large (3 simulations. All the fits have a good quality 
and the intercepts are zero within errors, with uncertainties of the order of 10~ 4 . 

3 Definition of the flow coupling 

Using our continuum result 



we define the gradient flow coupling for non-abelian gauge theories in the SF by means of 



This definition of the coupling is valid if the gauge field is coupled to fermions in arbitrary 
representations. As the reader may have noticed the scheme that defines the coupling 
depends not only on the quantities c,p,xq, but also on the value of the fermionic phase 
angle and the background field. In the simulations of the Schrodinger functional it is 
customary to include a phase angle 6 in the fermionic spatial boundary conditions. In 
principle different values of 9 are different schemes, although we have observed in some 
practical situations that the difference of the gradient flow coupling between 9 = and 
9 = 0.5 is below the 2%. 

Up to now we have worked exclusively with zero background fields, but the generaliza- 
tion to other values is straightforward. It only requires the modification of the heat kernels 
to preserve the value of the boundary fields and a modified form of the propagator [32]. 
Nevertheless common wisdom suggests that cutoff effects are reduced for zero background 
field, therefore we prefer to work in this scheme. In this case the definition of the coupling 
is also symmetric about xq = T/2 and we choose that value to minimize boundary effects. 
Also choosing p = T/L = 1 seems reasonable and leaves us with a one-parameter family of 
couplings, parametrized by the smoothing ratio c. 

By comparing the lattice and continuum behavior of the energy density as a function 
of c we can compute the leading order size of cutoff effects in the gradient flow coupling. 
As the reader can see in figure 2, the cutoff effects are large for small values of c, reach a 
minimum around c ~ 0.5 and then grow again. We recall that with c = 0.5 the smoothing 
radius is equal to L/2, and therefore one is effectively smoothing over all the lattice. For 



Af(c,p,x /T) 




(3.1) 



f GF (L) = [Af-\c, p, x /T) ■ t 2 (E(t, x Q ))] 



t=c 2 L 2 /S ■ 
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[hf-\c, p, x , a/L) - Af-^c, p, x )] /M~\c, p, x ) 



0.05 - 




-0.05 



-0.1 



0.3 0.35 0.4 0.45 0.5 

c 



Figure 2: Leading order relative cutoff effects of the gradient flow coupling as a function 
of the smoothing ratio c for different L/a at p = 1 and xq = T/2. 



c = 0.3 cutoff effects are smaller than 10% for a lattice of size L/a > 8, while for c = 0.4 
even the L/a = 6 lattice has cutoff effects of about 10%. 

This figure suggests using c = 0.5 as a preferred scheme, but later, when lattice simu- 
lations enter into the game, we will see that the statistical errors of the coupling also grows 
with c, and therefore in practice it is better to stay with c £ [0.3, 0.5], probably depending 
on the particular case, but this is the subject of the next section. 

We would also like to comment that if one is performing numerical simulations with 
the Wilson gauge action, one can benefit from smaller cutoff effects by using the lattice 
prediction to normalize the coupling. Defining 



fiV 2 -l)c 4 
M(c,p,x /T,a/L) = { — )C z 



128p 



P,P0 



p 2 cos 2 (api/2) - {pi cos(api/2)) 2 „ 2 



P 2 +Po 



+ 



p 2 cos 2 (ap /4) + \pl cos 2 (apj/2) 



P 2 +Pl 



Cp ( x ~ a / 2 ) 



the coupling would be given by 

g 2 GF (L) = \N-\c,p,x /T,a/L)-t 2 (E(t,x )) 



t=c 2 L 2 /8 



(3.3) 



(3.4) 



Obviously both definitions of the coupling differ only by cutoff effects. 

We finally want to mention that it is possible to define analogous couplings by using 
only the spatial components {Gik{t)Gik{t)) . In a lattice simulation one stays further away 
from the boundaries by not including plaquettes with links in the time direction. This may 
result in smaller cutoff effects, although this point needs further investigations. 
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0.1 0.2 0.3 0.4 0.5 0.6 0.3 0.35 0.4 0.45 0.5 

c c 

Figure 3: The gradient flow coupling as function of the flow time through c(t) = y/St/L 
for our lattices L\/a = 6, 8, 10, 12, 16 defined by a line of constant physics as described in 
the text. The right plot is a zoom of the dashed box in the left plot. The uncertainties are 
barely visible at this scale. 



IO) 4 




0.01 



c = 0.375 
c = 0.35625 
c = 0.3375 
c = 0.31875 

3 c = 0.3 
c = 0.28125 
c = 0.2625 
= 0.24375 
= 0.225 
c = 0.20625 
c = 0.1875 

c = 0.16875 
_l L 



_L_ 



0.02 

(a/Lf 



0.03 



0.04 




0.04 



Figure 4: A global view on the gradient flow coupling results for all five ensembles and 
c > c m ; n = max({a/L}). The connecting lines are drawn to guide the eye along results at 
constant smoothing ratio with values given in the plot. 



4 Non-perturbative tests 

In this section we would like to analyze the gradient flow coupling numerically. We want 
to estimate both the size of cutoff effects and the numerical cost of evaluating the new 
gradient flow coupling. The main result of this section is that both quantities depend on 
the particular scheme via the parameter c. When c is increased cutoff effects decrease, but 
the numerical cost increases. We find that the window of values c € [0.3, 0.5] allows a very 
precise determination with a mild continuum extrapolation. 

4.1 Line of constant physics 

As framework for our tests we choose a set of iVf = 2 Schrodinger functional simulations 
at a line of constant physics as given through 

ffSF^i) = u = 4 - 484 and m ( L i) = i 
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where ^g F is the traditional SF coupling and m the renormalized PCAC mass. From 
reference [34] we know that the physical volume is roughly L\ ~ 0.4 fm. The available 
five different ensembles are lattices with L/a £ {6, 8, 10, 12, 16} at T = L with vanishing 
boundary gauge fields and a fermionic phase angle 6 = 0.5. Each ensemble consists of at 
least 8000 configurations separated by r meas = 10 molecular dynamic units (MDU). We 
refer the reader to the appendices in [35] for any unexplained detail concerning the physics 
and run parameters. 

We would like to measure the value of the gradient flow coupling in these ensembles. 
Since they have been tuned to have constant SF coupling, equivalent to constant volume 
in the continuum, we define the function 

r ~ , o 1 «=4.484, m sca =0, 9=0.5 

n(u; c, a/L) = M~\c, 1, 1/2, a/L) ■ t 2 (E(t, T/2)) (4.2) 

L J t=c i L i /8 

that at fixed c has the gradient flow coupling as continuum limit 

5q F (L) = cj(u; c) = lim Q(u;c,a/L). (4-3) 
a/L— s-o 

We will also use some ensembles with larger lattices {L\/a = 24,32,40) but lower 
statistics. These have been defined by a slightly different line of constant physics: 

5§ F (Li/2) = u = 2.989 and m(Li/2) = , (4.4) 

that is, at a fixed SF coupling corresponding to half the scale L\. Both LCP's are related 
through the step scaling function in two-flavour QCD [8], 

g\-p{L x ) = cr(2.989) = 4.484(48) , (4.5) 

and thus differ only by cutoff effects. The statistics for this second set of ensembles is 
smaller (~ 800 measurements). For these additional lattices we define a function similar 
to O according to 

~ r * 1 r, 1 u=a(u), m SGa =0, 0=0.5 

n(«; c, a/L) = Af- l (c, 1,1/2, a/L) ■ t 2 (E(t,T/2)) , (4.6) 

L J t=c^L / /8 

that also has the same continuum limit. 

All ensembles have been tuned to have constant SF coupling only with some statis- 
tical accuracy. This propagates into an uncertainty in the determination of the functions 
Cl(u; c, a/L) and ft(u; c, a/L). This error can be estimated by simple propagation of errors, 
for example applied to Q(u;c,a/L) we have 

an 



SU(u; c, a/L) 



du 



8u. (4.7) 



To evaluate this uncertainty we use another ensemble (labeled 12* in table 2) with a 
slightly different value of (3 but also tuned to have vanishing quark mass. By evaluating 
both the SF coupling and O on this ensemble we can numerically estimate the derivative 
in equation (4.7) 2 . This source of error in fact dominates the error budget of Q(u;c,a/L), 
which anticipates that the new coupling is numerically more precise. 



2 For c = 0.3, 0.4, 0.5 we obtain (dQ/du) = 0.7, 1.1, 1.5 respectively. 
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L/a 


6 


8 


10 


12 


16 


12* 


P 


5.2638 


5.4689 


5.6190 


5.7580 


5.9631 


5.8120 


^sca 


0.135985 


0.136700 


0.136785 


0.136623 


0.136422 


0.136617 


mcas 


12160 


8320 


8192 


8280 


8460 


2392 




4.423(75) 


4.473(83) 


4.49(10) 


4.501(91) 


4.40(10) 


4.218(49) 


Q(u: 0.3, a/L) 


4.818(55) 


4.728(61) 


4.627(73) 


4.518(67) 


4.441(74) 




9gf 


4.8178(46) 


4.7278(46) 


4.6269(47) 


4.5176(47) 


4.4410(53) 


4.310(8) 


Tint 


0.57(2) 


0.51(2) 


0.62(3) 


0.66(3) 


0.92(6) 


0.67(6) 


Rns x 10 3 


0.95(2) 


0.97(2) 


1.02(3) 


1.04(3) 


1.20(4) 


- 


Sl(u;OA,a/L) 


6.009(80) 


5.699(88) 


5.60(11) 


5.484(96) 


5.41(11) 




9 

(J™ 


6.0090(86) 


5.6985(86) 


5.5976(97) 


5.4837(97) 


5.410(12) 


5.182(16) 


Tint 


0.55(2) 


0.52(2) 


0.70(4) 


0.76(4) 


1.24(9) 


0.73(7) 


i?NS X 10 3 


1.43(3) 


1.51(4) 


1.73(5) 


1.77(5) 


2.23(9) 




n(u;0.5,o/L) 


7.11(12) 


6.82(13) 


6.76(15) 


6.66(14) 


6.60(15) 






7.106(14) 


6.817(15) 


6.761(19) 


6.658(19) 


6.602(24) 


6.223(29) 


Tint 


0.54(2) 


0.57(2) 


0.82(5) 


0.89(5) 


1.49(12) 


0.82(9) 


i?NS x 10 3 


1.97(4) 


2.26(5) 


2.85(9) 


2.81(9) 


3.6(2) 





Table 2: Lattice run parameters (see [35] for more details) and results for the gradient 
flow coupling at smoothing ratio c G {0.3,0.4,0.5}. Errors are computed using the T- 
method [36]. We show the values of the gradient flow coupling g|, F (Li) and of the function 
Cl(u;c,a/L). Furthermore, we quote the integrated auto-correlation time T mt of g|, F (Li), 
estimated in units of the measurement frequency r mcas = 10 MDU (which is the same for 
each lattice), and the noise-to-signal ratio Rns- The lattice labeled 12* is used to estimate 
(dU/du). 



4.2 Numerical results and computing cost 

Figure 3 shows the gradient flow coupling as a function of c for the different ensembles. 
For c G [0.3,0.5] we observe a monotonous behavior of ^|. F with a/L. As figure 4 shows 
this seems to be the scaling region of the gradient flow coupling for lattices with L/a > 8, 
and therefore this is the region in which we will focus from now on. 
In figure 5 we present the noise-to-signal ratio 

i?NS = ^ (4.8) 

as obtained in our analysis as function of c for the individual lattices. We observe that the 
noise-to-signal ratio increases with increasing c (see table 2). Although our statistics does 
not allow us to draw definite conclusions, we observe a behavior compatible with a power- 
like scaling of i?NS with c. The behavior seems to be universal and independent of a/L in 
contrast to the traditional SF coupling, that has a divergent variance when approaching 
the continuum [37]. The product V^mcas-^NS can directly be translated in the cost of 
obtaining the new gradient flow coupling with some precision. 8000 measurements are 
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L/a 


24 


32 


40 


a 
P 


a O/iQQ 
0.Z4OO 


0.40 / 4 


a coon 
O.O00U 


^sca 


fl 1 ^^Ql D/l 
U. loOy 1U4 


U. loOOZlU 


fi iqci 090 

u.iooiyzo 


-^mcas 


632 


800 


850 


Tineas 


10 MDU 


10 MDU 


4 MDU 


7J2 IT lcs\ 








n(u;0.3,a/L) 


4.405(45) 


4.402(53) 


4.335(71) 


772 

9gf 




A AC\0(A(X\ 




7*int / ^"meas 


9 nf^ 


3 on 9"i 


1 9C e i , \ 


Q(u;OA,a/L) 


5.39(8) 


5.46(12) 


5.34(16) 


9gf 


5.39(6) 


5.46(11) 


5.34(15) 


7"int 1 tineas 


2.5(7) 


6.1(2.1) 


19(9) 


Q(u; 0.5, a/L) 


6.64(14) 


6.87(25) 


6.67(30) 


9gf 


6.62(12) 


6.87(24) 


6.67(29) 


T\jxt 1 7"meas 


2.9(9) 


7.4(2.8) 


21(10) 



Table 3: Same as table 2 but for the second line of constant physics, eq. (4.4). 
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Figure 5: Numerical cost of evaluating the gradient flow coupling. In the plot we show 
V NmeasRNS as a function of c. For the interesting values of c > 0.3 there seem to be a 
power law scaling roughly universal for all values of L/a. 



enough to achieve a precision of 0.1% for c = 0.3, while for c = 0.5 the precision decreases 
to 0.35%. 

In figure 6 we plot the integrated auto-correlation time Tint f° r the different ensembles. 
The lattices L/a = 6,8 have T; n t ~ 0.5 which means that the available configurations 
are too far separated in Monte Carlo simulation time to detect any auto-correlations in 
the chain. The L/a = 10, 12, 16 lattices show a clear increase of auto-correlations with 
increasing L/a. As can be inferred from figure 6 this increase is compatible with a scaling 
~ l/a 2 towards the continuum. This is much in accordance with the conjecture of [38, 39] 
for the scaling behaviour of the HMC algorithm in an interacting theory, since we do not 
expect non-zero topological charge sectors to contribute significantly at this small physical 
volume. 
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Figure 6: (Left) Integrated auto-correlation time as a function of c for the lattices L/a 
6,8,10,12,16. (Right) (a/ L) 2 T in t as a function of L/a for three representative cases c 
0.3,0.4,0.5 and all lattices up to L/a = 40. 



4.3 Cutoff effects 

For the continuum approach of the function Q(u;c,a/L) (and also of Q(u;c,a/L)) we 
observe a behaviour dominated by a linear scaling in (a/L) 2 . Hence, we choose 

n(u;c,a/L) = lo(u; c) { 1 + A(u; c) ■ (a/L) 2 } (4.9) 

as fit ansatz to extract the continuum limit u(u; c) = <7q F (L) and the leading cutoff effects 
A(u;c). Figure 7 shows examples of these extrapolations for three representative cases 
c = {0.3, 0.4, 0.5}. We observe that cutoff effects decrease with increasing c. Quantitatively 
this can be estimated by looking at the relative size of cutoff effects through the ratio 

R(u; c; a/L) = 11 , 4.10 

uj(u; c) 

that decreases by a factor 2 when c is increased from 0.3 to 0.5 (see figure 7). 

Since in general cutoff effects of the SF step scaling function are very small, and given 
the fact that cutoff effects of £l(u;c,a/L) change with c, we think that the cutoff effects 
in Q(u;c, a/L) and f2(u;c, a/L) are dominated by the lattice spacing dependence of the 
new gradient flow coupling. Therefore we expect both functions to show roughly the same 
scaling behavior. Although the points corresponding to f2(u; c,a/L) have not been used for 
the previous continuum extrapolations, we have added the points to the plots in figure 7. 
The data fits well into the expected scaling behavior. 

Summarizing figure 7 and table 2 we can say that when c is increased from c = 0.3 
to c = 0.5 the relative cutoff effects decrease by about a factor of 2. We remember from 
previous sections that this decrease of cutoff effects comes at the expense of an increased 
relative statistical error (about three times larger when going from c = 0.3 to c = 0.5). 

4.4 Mass dependence 

Last but not least we have studied the mass dependence of the gradient flow coupling 7/q F as 
it was done for ^g F in [35]. For this purpose we have generated an ensemble with L/a = 8 
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Figure 7: (Left) Continuum extrapolations of Q(u; c, a/L) for three different values of c = 
0.3, 0.4, 0.5. Only the lattices with L/a = 8, 10, 12, 16 are used for the fit, but the L/a = 6 
lattice as well as the larger lattices L/a = 24, 32, 40 are in the plot. (Top right) The value 
of the gradient flow coupling in the continuum as a function of c. (Bottom right) Relative 
size of the continuum extrapolation for the three representative cases c = 0.3, 0.4, 0.5. 



at the same value of the bare coupling as the available one, but with a non-zero quark 
mass. Actually the bare parameters of the simulation correspond to the lattice labeled as 
8* in [35], and the interested reader is encouraged to consult the original work for more 
details. Defining the dimensionless PCAC quark mass z = Lm, we obtain 



^GF 



dz 



u =4.484 



0.19(7) for c = 0.3 
0.17(9) for c = 0.4 



(4.11) 



to be compared to the corresponding value of 1.4(4) for the Schrodinger functional coupling. 
The mass dependence of the gradient flow coupling as defined in the present paper is smaller 
by an order of magnitude. 



5 Conclusions 

The gradient flow can be used to defined a renormalized coupling at a scale fi = In 
this work we have studied the perturbative behavior of the gradient flow in the Schrodinger 
functional. By setting the renormalization scale proportional to the linear size of the SF 
box, n = 1 /\^8t = 1 /cL, we have defined a family of running coupling constants valid for an 
arbitrary SU(N) gauge field coupled to arbitrary fermions. Since this coupling definition 



- 16 - 



does not depend on any scale but the finite volume, it can be used for finite-size scaling 
purposes. 

The coupling constant can be defined for different values of the background field in 
the SF. Since one expects cutoff effects to be smaller for the case of vanishing background 
field this is the case that we have studied in more detail. From our perturbative analysis 
we have been able to study the size of cutoff effects to leading order. Cutoff effects tend 
to be relatively large for either small or large values of c, but very mild for c £ [0.3,0.5]. 
As an example for c = 0.3 the difference between the coupling in a L/a = 8 lattice and the 
continuum is around 10%, while for c = 0.5 the difference between the continuum and the 
value in an L/a = 6 lattice is below 4%. 

We have analyzed a total of five ensembles tuned to have a constant SF coupling in a 
physical volume of L ~ 0.4 fm. The cutoff effects observed in the gradient flow coupling 
on these ensembles shows a similar overall behavior as expected from earlier perturbative 
considerations. Therefore, we conclude that the gradient flow coupling has mild cutoff 
effects if c is chosen wisely even for small lattices. 

We have analyzed the numerical cost of evaluating the gradient flow coupling, and 
see that it increases with c. Nevertheless it can be computed very precisely with a modest 
numerical effort. For c = 0.3 we find that roughly 8000 measurements are enough to obtain 
a precision of 0.1% even in our larger lattice L/a = 16, while for the worst analyzed case 
at c = 0.5 this precision drops to the still very good figure of about 0.35%. 

What is the most convenient scheme? From a practical point of view there is no need 
to settle this discussion - and thus the unique definition of the scheme - immediately since 
one measures the gradient flow coupling at different values of c in any case while integrating 
the Wilson flow. Probably the optimal choice will depend on the particular problem. 

In summary we conclude that the gradient flow coupling in the SF has several prac- 
tical advantages. It is valid for arbitrary SU(N) gauge fields coupled to fermions in any 
representation, having the universal two loop beta function. It has mild cutoff effects and 
can be cheaply evaluated numerically with high precision. It has a very small dependence 
on quark masses and is naturally defined with vanishing background field, avoiding the 
generation of new ensembles just to compute the running coupling. 

We think that the computation of the next order in the perturbative behavior of the 
gradient flow coupling in the SF is very interesting, in particular in connection to the 
determination of the Lambda parameter. This computation can also shed some light on an 
optimal value of c for the matching between different schemes. There are many applications 
for this new coupling, but we are particularly interested in using it to compute the step 
scaling function of QCD. We believe that this can be achieved with a very high accuracy 
with a modest computational effort. 

Acknowledgments 

In particular, the authors would like to thank R. Sommer for suggesting us (with immense 
patience) to look into this very interesting problem and also for guiding us in the steps of 
this work. 



-17- 



We have profited from the generosity of many individuals who have provided us with 
code for several purposes. We thank M. Liischer for making available a code to integrate 
the flow equations numerically, to H. Simma for the invaluable help in adapting the code 
to our needs, and finally to M. Marinkovic (SF-HMC) and the openQCD team (M. Liischer, 
S. Schaefer and J. Bulava for implementing SF boundary conditions) for their very useful 
codes that we have used at different stages of this project to generate new ensembles. 

A. R. also wants to thank A. Gonzalez- Arroyo, S. Sint, U. Wolff and F. Virotta for the 
many illuminating discussions concerning the role of the boundary conditions, the SF and 
the Wilson flow. Thanks, guys! 

This work is supported by the Deutsche Forschungsgemeinschaft in the SFB/TR 09 
and the John von Neumann institute for computing. Our numerical computations were 
performed on the PAX cluster at DESY, Zeuthen. 

A Notation 

Momenta are always defined to be 



2im; 

Pi = -jf (A.l) 



for the periodic, spatial directions and 

27rn ° ( A 

Po = ~y~ ( A - 2 ) 
in the time direction (xq). In the continuum sum over momenta are abbreviated by 

• i 2-7rn 2irn 

= 2^ > Wlth p = ~ OT p = — ( A - 3 ) 

p n=—oo 

while in the lattice we have finite sums 

Lj1a-\ 

with T)„- = - 

L 



^2= ' with = ' ( A - 4 ) 

Pi rii=—L/2a 

£= ^ , withp = ^. (A.5) 

PO n =-T/a 

If we impose Dirichlet/Neumann boundary conditions on the interval xq £ [0,T] 

/(seo)Uo=o,t = Dirichlet, (A. 6) 

9x0/(20)1x0=0^ = Neumann, (A.7) 

the Laplacian has eigenfunctions and eigenvalues given by 

f 7Tn o x o \ 

Sp { x o) = sm(p x /2) = sin y j , (A.8) 

Cp (xQ) = cos(p x /2) = cos I J , (A.9) 
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respectively. The corresponding eigenvalues are given through 

2 



d l s po( x o) = ~ (y) s po( x o), 

dl c P0 {x ) = - (y) c po (x ) , 
and these functions obey the completeness relations 

s p( x ) s p( x ) = ^xx' i 

v 

y ] c p( x ) c p( x ) = Sxx' i 

P 

1 f T 1 

2T J dx s p( x ) s i( x ) = 2 ^ 



1 f T , 1 



dxc p (x)c g (x) = ^ (S pq + 5 P . 

On the lattice derivatives are substituted by finite differences 

d x f(x) = -(f( x + a )-f(x)) , 
a 

d*f(x) = -(/(*)-/(* -a)) , 
a 

and defining the family of lattice momenta 

V i i = - sin I a^f J , = - sin (ajp M ) , p M = - sin I a-f 

a V 2 / a a V 4 

the discrete Laplacian d x d* has eigenvalues/eigenfunctions with periodic boundary 
tions 

8d*e ipx = -p 2 e ipx , 
where p = 2nn/L. With Dirichlet/Neumann boundary conditions one has 

dd*s p (x) = —p 2 s p (x) Dirichlet , 
dd*c p (x) = —p 2 dp(x) Neumann, 

where now 

'px 



Sp(x) = sin ^ 

„ f s \p(x + a/2)~ 
Cp{x) = cos 
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satisfy similar completeness relations 

v 

1 r ~ 1 

— s p (x)s q (x) 



x=0 



2 (3pq 3p,—q) 



v 

1 T_1 



(Spq + £p, 



(A.24) 

(A.25) 
(A.26) 

(A.27) 



x=0 



B Heat Kernels 



Here we will review some known properties of heat kernels. The interested reader will like 
to read appendix C of [40], with more information of heat kernels and the SF. 

B.l Continuum heat kernels 

Heat kernels are fundamental solutions (i.e. a solution that is a delta source at t = 0) to 
the heat equation, that in one dimension reads 



d d 2 



(B.l) 



To obtain a solution to a given problem one has to choose appropriate boundary conditions. 
In the case of Euclidean space the heat kernel is well known to be 

(x-x') 21 



K(x,x',t) = (47rt)- 1/2 exp 



4/ 



(B.2) 



Since the heat equation is a linear equation one can construct heat kernels with different 
boundary conditions with a linear combinations of the heat kernel in Euclidean space. In 
particular 

(x - x' + nL) 21 



it 



(B.3) 



K p (x,x',t) = (4vrt)- 1 / 2 ex P 

n=— oo 

is, by construction, both a solution to the heat equation and is periodic in x with period 
L. Therefore it corresponds to the heat kernel in S . Following a similar reasoning it is 
easy to see that 

(x - x' + 2nT) 2 " 



K D (x,x',t) = (4irt)~ 1/2 i ex P 



4/ 



exp 



(x + x 1 + 2nTY 
it 



(B.4) 



and 



K N (x,x',t) = (4nt)^ 2 f; {ex P ( >-*' + 2raT)2 )+ex P (- 



(x + x' + 2nTf 



(B.5) 
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are heat kernels with Dirichlet and Neumann boundary conditions on the interval [0, T] 
respectively. It is worth mentioning that, via the Jacobi imaginary transformation, the 
periodic heat kernel is nothing more than the third Jacobi theta function 3 

K p (x,x',t) = itf 3 (l(x-J)\*p) = \ ± e-mVPO-O. (B.6) 

^ ' n=— oo 

This last expression, that can also be obtained using Poisson summation formula, turns 
out to be convenient for our computations. 

B.2 Discrete heat kernels 

When performing computations of the Wilson flow on the lattice we find a type of heat 
equation in which the time variable is continuous but the Laplacian is substituted by a 
discrete version 

d t f(x,t) = d x B*J(x,t) (B.7) 

Taking appropriate boundary conditions into account, we call fundamental solutions of this 
equation discrete heat kernels. The most easy way to construct them is by noting that one 
can formally write 

K = exp\td x d*} (B.8) 

and now the task of finding the heat kernels is reduced to finding eigenvalues /eigenfunctions 
of the discrete Laplacian with the correct boundary conditions. By recalling the notions of 
appendix A one can immediately write 



(x,x',t) = I^ e -P 2 V^-*'), (B.9) 
v 

K D (x,x',t) = I^ e ^ 2 %(x)5 p (:r'), (B.10) 
p 

K N (x, x', t) = ^J2 ^ p \(x)c p (x') . (B.ll) 



B.3 Properties 



The following properties are straightforward, but fundamental to easily reproduce our 
results: 

/ dx' K p (x,x',t)e ipx ' = e~ p2t e ipx , (B.12) 
J o 

f dx' K D (x,x',t)s p (x') = e-(i) 2i Sp(x) , (B.13) 
Jo 

f T A r dK D { X ,x',t) , _(E) 2 t , s rT3 iA, 

/ dx — —cJx ) = zpe t c p (x) , (B.14) 

Jo dx 

dx' K N (x, x', t)cp(x') = e-(§)\p(x) . (B.15) 



3 A standard reference with the same conventions used here is [41] 
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These relations have discrete analogons 

L-l 



K p (x, x', t)e ipx ' = e* V px , (B.16) 

x'=0 
T 

K D (x, x', t)s p (x') = e- p2t s p (x) , (B.17) 

x'=0 

y dK D (x,x',t) , = ^ cos(ap/2)e -p 2 t c , x s m lg) 
ox 

J2 K N {x, x', t)c p {x') = e- pH Cp(x) . (B.19) 



x'=0 
T 



x'=Q 

Finally we note that the construction of heat kernels in more than one dimension is done 
simply by multiplying one-dimensional heat kernels. For example 

3 

K SF {x,x',t) =l[k p (x i ,x' i ,t)K D (x ,x' ,t) (B.20) 

i=i 

is a 4-dimensional heat kernel with periodic b.c. in the three space dimensions and Dirichlet 
b.c. in the time dimension. 

C Propagators 

In the SF the free gluon propagator on the lattice is defined as 

(A a ^p,x )A b l/ {q,y )} = L 3 £ a6 5 Pi _ q ZV(p,x ,yo) ( ai ) 
and in the Feynman gauge given by 

D ik (p,x ,y ) = 5 ik d(p,x ,y ) , (C.2) 
D ok (p, x , y ) = D k0 (p, x ,y ) = , (C.3) 
Doo(p, x , y ) = n(p, x ,y ) . (C.4) 

The functions d(p, XQ,yo) and n(p, XQ,yo) are defined as 

d p, x , y = - > - 2 . 2 — , (C.5 



Po 

n(p,x ,y = - > (P^O). C.6 

T^f P 2 +P 2 



The reader interested in the explicit form of the propagator for different background 
fields, or in different gauges, should consult the original references [31, 32] 

To obtain the propagator in the continuum it is enough to take the continuum limit 
of the previous expressions, that simply amounts to taking out the hats, and making the 
sums run trough all integers. 
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D Adaptive size integrators for the Wilson flow 



On the lattice the flow equation has the form 

Z(V)VJx,t). (D.l) 



2 dV fl (x,t) 
1 dt ~ 



Following the advice in [27], we use the third order Runge-Kutta scheme given by 

W Q = V lt (x,t), (D.2) 

W x = e X p|iz |wo, (D.3) 

W 2 = exp |^Zi - ^Z | Wt , (D.4) 

V^x,t + a 2 e) =exp|^ 2 -^Z 1 + ^Z |w 2 , (D.5) 



where 



Zi = eZ{W t ) . (D.6) 



We simply want to point out that with one extra computation, one can get a second 
estimate of V /1 (x, t + a?e) 

V^x, t + a 2 e) = exp {-Z + 2Z X } W 2 . (D.7) 

This last scheme is of second order, as the reader can easily check. Given two N x N 
complex matrices A, B, one can define a distance between them as follows 

This distance can be used to estimate the error made by the lower order integrator 

d = max {distO^x.t + a 2 e),V^(x,t + a 2 e))} . (D.9) 

With this information one can adjust the step size so that the error in the integration 
does not exceed certain tolerance 5. After each integration step e is updated according to 

e— $-60.95^! (D.10) 

and if d > 5 the integration step is repeated. 

The reason why this scheme is efficient for the particular case of the Wilson flow is 
related to its smoothing properties. Close to t = the configuration is rough and therefore 
a very fine integration is needed. But as t increases the configuration is more and more 
smooth, and one can have a very precise integration with a large e. In practical cases it has 
saved us around a factor 4 in the number of integration steps. It also has the advantage 
that e is tunned automatically, and one does not need to worry about the step size, but 
only plug in the desired tolerance. 
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